Hands-on -Exercise 4

Spatial Point Pattern Analysis (SPPA) is the evaluation of the pattern or distribution of a set of points on a surface. The points may represent:

events such as crimes, traffic accidents, or disease onsets, or business services (e.g., coffee shops and fast-food outlets) or facilities such as childcare centres and eldercare centres. First-order Spatial Point Pattern Analysis (1st-SPPA) focuses on understanding the intensity or density of points across a study area. It examines how the distribution of points varies over space, essentially identifying trends or patterns in point density. This type of analysis deals with the individual locations of points and their distribution, without considering interactions between them.

In essence, 1st-SPPA helps answer questions such as:

Where are points most densely located within the study area? Is point density uniform, or does it vary across space? How spread out is the point pattern? In this chapter, you will gain hands-on experience with spatstat to perform two commonly used 1st-SPPA methods:

The specific questions we would like to answer are as follows:

Are the childcare centres in Singapore randomly distributed throughout the country? If the answer is not, then the next logical question is where are the locations with higher concentration of childcare centres?

Installing and Loading the R packages

pacman::p_load(sf, terra, spatstat, 
               tmap, rvest, tidyverse)

Importing and Wrangling Geospatial Data Sets

mpsz_sf <- sf::st_read("../Hands_on_Ex02/Data/geospatial/MasterPlan2019SubzoneBoundaryNoSeaKML.kml") |>
  sf::st_zm(drop = TRUE, what = "ZM") |>
  sf::st_transform(crs = 3414)

Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source 
  `C:\Users\kchan\Desktop\chandru-ko\ISSS608-VAA\Hands-on_Ex\Hands_on_Ex02\Data\geospatial\MasterPlan2019SubzoneBoundaryNoSeaKML.kml' 
  using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

Next, build a function called extract_kml_field for extracting REGION_N, PLN_AREA_N, SUBZONE_N, SUBZONE_C from Description field by using the code chunk below.

extract_kml_field <- function(html_text, field_name) {
  if (is.na(html_text) || html_text == "") return(NA_character_)
  
  page <- read_html(html_text)
  rows <- page %>% html_elements("tr")
  
  value <- rows %>%
    keep(~ html_text2(html_element(.x, "th")) == field_name) %>%
    html_element("td") %>%
    html_text2()
  
  if (length(value) == 0) NA_character_ else value
}

mpsz_sf <- mpsz_sf %>%
  mutate(
    REGION_N = map_chr(Description, extract_kml_field, "REGION_N"),
    PLN_AREA_N = map_chr(Description, extract_kml_field, "PLN_AREA_N"),
    SUBZONE_N = map_chr(Description, extract_kml_field, "SUBZONE_N"),
    SUBZONE_C = map_chr(Description, extract_kml_field, "SUBZONE_C")
  ) %>%
  select(-Name, -Description) %>%
  relocate(geometry, .after = last_col())

mpsz_cl <- mpsz_sf %>%
  filter(SUBZONE_N != "SOUTHERN GROUP",
         PLN_AREA_N != "WESTERN ISLANDS",
         PLN_AREA_N != "NORTH-EASTERN ISLANDS")

write_rds(mpsz_cl, 
          "Data/mpsz_cl.rds")

Next, code chunk below will be used to import the childcare Service data downloaded from data.gov.sg into R environment as sf data frame called chilcare_sf. The Simple Feature Geometry (sfg) of this geospatial data layer if

childcare_sf <- st_read("Data/ChildCareServices.kml") %>% 
  st_zm(drop = TRUE, what = "ZM") %>%
  st_transform(crs = 3414)

Reading layer `CHILDCARE' from data source 
  `C:\Users\kchan\Desktop\chandru-ko\ISSS608-VAA\Hands-on_Ex\Hands_on_Ex04\Data\ChildCareServices.kml' 
  using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

Mapping the geospatial data sets

After checking the referencing system of each geospatial data data frame, it is also useful for us to plot a map to show their spatial patterns.

tmap_mode('view')

ℹ tmap mode set to "view".

tm_shape(childcare_sf)+
  tm_dots()

Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite

tmap_mode('plot')

ℹ tmap mode set to "plot".

Geospatial Data wrangling

Converting sf data frames to ppp class

spatstat requires the point event data in ppp object form. The code chunk below uses [as.ppp()] of spatstat package to convert childcare_sf to ppp format.

childcare_ppp <- as.ppp(childcare_sf)

Next, class() of Base R will be used to verify the object class of childcare_ppp.

class(childcare_ppp)

[1] "ppp"

You can take a quick look at the summary statistics of the newly converted ppp object by using the code chunk below.

summary(childcare_ppp)

Marked planar point pattern:  1925 points
Average intensity 2.417323e-06 points per square unit

Coordinates are given to 11 decimal places

Mark variables: Name, Description
Summary:
     Name           Description       
 Length:1925        Length:1925       
 Class :character   Class :character  
 Mode  :character   Mode  :character  

Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
                    (33590 x 23700 units)
Window area = 796335000 square units

Creating owin object

When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.

The code chunk below, as.owin() of spatstat is used to covert mpsz_sf into owin object of spatstat.

sg_owin <- as.owin(mpsz_cl)

Again, class() will be used to verify the object class of sg_owin.

class(sg_owin)

[1] "owin"

The result above confirmed that sg_owin is indeed in owin object.

sg_owin object can be displayed by using plot() function.

plot(sg_owin)

Combining point events object and owin object

In this last step of geospatial data wrangling, we will extract childcare events that are located within Singapore by using the code chunk below.

childcareSG_ppp = childcare_ppp[sg_owin]

Nearest Neighbor Analysis (NNA) is a spatial statistics method that calculates the average distance between each point and its closest neighbor to determine if a pattern of points is clustered, dispersed, or randomly distributed.

Clark-Evans test is a specific statistical method used within NNA to quantify whether a point pattern is clustered, random, or uniformly spaced, using the Clark-Evans aggregation index (R) to describe this pattern. NNA provides a numerical value that describes the degree of clustering or regularity, and the Clark-Evans test calculates a specific index (R) for this purpose

Perform the Clark-Evans test without CSR

clarkevans.test() of spatstat.explore package support two Clark-Evans test, namely: without CRS and with CRS. In the code chunk below, Clark-Evans test without CSR method is used.

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"))

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.53532, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

Perform the Clark-Evans test with CSR

In the code chunk below, the argument method = “MonteCarlo” is used. In this case, the p-value for the test is computed by comparing the observed value of R to the results obtained from nsim (i.e. 39, 99, 999) simulated realisations of Complete Spatial Randomness conditional on the observed number of points.

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                method="MonteCarlo",
                nsim=99)

    Clark-Evans test
    No edge correction
    Monte Carlo test based on 99 simulations of CSR with fixed n

data:  childcareSG_ppp
R = 0.53532, p-value = 0.01
alternative hypothesis: clustered (R < 1)

Kernel Density Estimation Method

Kernel Density Estimation (KDE) is a valuable tool for visualising and analyzing first-order spatial point patterns. It is widely considered a method within Exploratory Spatial Data Analysis (ESDA) because it’s used to visualize and understand spatial data patterns by transforms discrete point data (like locations of childcare service, crime incidents or disease cases) into continuous density surfaces that reveal clusters and variations in event occurrences, without making prior assumptions about data distribution. It helps to begin understanding data distribution, identify hotspots, and explore relationships between spatial variables before performing more rigorous analysis.

Working with automatic bandwidth selection method The code chunk below computes a kernel density by using the following configurations of density() of spatstat:

bw.diggle() automatic bandwidth selection method. Other recommended methods are bw.CvL(), bw.scott() or bw.ppl(). The smoothing kernel used is gaussian, which is the default. Other smoothing methods are: “epanechnikov”, “quartic” or “disc”. The intensity estimate is corrected for edge effect bias by using method described by Jones (1993) and Diggle (2010, equation 18.9). The default is FALSE.

kde_SG_diggle <- density(
  childcareSG_ppp,
  sigma=bw.diggle,
  edge=TRUE,
  kernel="gaussian") 

The plot() function of Base R is then used to display the kernel density derived.

plot(kde_SG_diggle)

In the code chunk below, summary() of Base R is used to print the summary report of the

summary(kde_SG_diggle)

real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [2667.538, 55941.94] x [21448.47, 50256.33] units
dimensions of each pixel: 416 x 225.0614 units
Image is defined on a subset of the rectangular grid
Subset area = 669941961.12249 square units
Subset area fraction = 0.437
Pixel values (inside window):
    range = [-6.584123e-21, 3.063698e-05]
    integral = 1927.788
    mean = 2.877545e-06

Before we move on to next section, it is good to know that you can retrieve the bandwidth used to compute the kde layer by using the code chunk below.

bw <- bw.diggle(childcareSG_ppp)
bw

   sigma 
295.9712 

Rescalling KDE values

In the code chunk below, rescale.ppp() is used to covert the unit of measurement from meter to kilometer.

childcareSG_ppp_km <- rescale.ppp(
  childcareSG_ppp, 1000, "km")

Now, we can re-run density() using the resale data set and plot the output kde map.

kde_childcareSG_km <- density(childcareSG_ppp_km,
                              sigma=bw.diggle,
                              edge=TRUE,
                              kernel="gaussian")

Next, plot() is used to plot the kde object as shown below.

plot(kde_childcareSG_km)

Working with different automatic badwidth methods

Beside bw.diggle(), there are three other spatstat functions can be used to determine the bandwidth, they are: bw.CvL(), bw.scott(), and bw.ppl().

Let us take a look at the bandwidth return by these automatic bandwidth calculation methods by using the code chunk below.

bw.CvL(childcareSG_ppp_km)

   sigma 
4.357209 

bw.scott(childcareSG_ppp_km)

 sigma.x  sigma.y 
2.159749 1.396455 

bw.ppl(childcareSG_ppp_km)

   sigma 
0.378997 

bw.diggle(childcareSG_ppp_km)

    sigma 
0.2959712 

Baddeley et. (2016) suggested the use of the bw.ppl() algorithm because past experience shown that it tends to produce the more appropriate values when the pattern consists predominantly of tight clusters. But they also insist that if the purpose of once study is to detect a single tight cluster in the midst of random noise then the bw.diggle() method seems to work best.

The code chunk beow will be used to compare the output of using bw.diggle and bw.ppl methods.

kde_childcareSG.ppl <- density(childcareSG_ppp_km, 
                               sigma=bw.ppl, 
                               edge=TRUE,
                               kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG_km, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")

Working with different kernel methods

By default, the kernel method used in density.ppp() is gaussian. But there are three other options, namely: Epanechnikov, Quartic and Dics.

The code chunk below will be used to compute three more kernel density estimations by using these three kernel function.

par(mfrow=c(2,2))
plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")
plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")
plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")
plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")

Fixed and Adaptive KDE

Computing KDE by using fixed bandwidth

Next, you will compute a KDE layer by defining a bandwidth of 600 meter. Notice that in the code chunk below, the sigma value used is 0.6. This is because the unit of measurement of childcareSG_ppp_km object is in kilometer, hence the 600m is 0.6km.

kde_childcareSG_fb <- density(childcareSG_ppp_km,
                              sigma=0.6, 
                              edge=TRUE,
                              kernel="gaussian")
plot(kde_childcareSG_fb)

Computing KDE by using adaptive bandwidth

Fixed bandwidth method is very sensitive to highly skew distribution of spatial point patterns over geographical units for example urban versus rural. One way to overcome this problem is by using adaptive bandwidth instead.

In this section, you will learn how to derive adaptive kernel density estimation by using density.adaptive() of spatstat.

kde_childcareSG_ab <- adaptive.density(
  childcareSG_ppp_km, 
  method="kernel")
plot(kde_childcareSG_ab)